## ----setseedch9, echo = FALSE--------------------------------------------
set.seed(1982)

knitr::opts_chunk$set(
                      echo = FALSE,
                      results = "hide",
                      warning = FALSE,
                      message = FALSE,
                      error = FALSE
)


## ----loadlibsch9---------------------------------------------------------
library(dplyr)
library(ggplot2)
library(reshape2)
library(zoo)
library(memisc)
library(pander)
library(pscl)
library(rstanarm)
library(stringr)
library(Amelia)
source("common_funcs.R")


## ----loaddatach9---------------------------------------------------------
load("ch2_data.RData")

## From the codebook:

## Description: This variable records the outcome (or outcomes) of the case. 

## Valid values: Valid values are 

## - A 
## - B
## - C
## - D
## - E

## corresponding to:

## - The appeal was allowed
## - The appeal was allowed in part and dismised in part
## - The appeal was dismissed
## - The reference questions were answered
## - Other order
## Valid values may be joined together with a semi-colon, but no space. "A;C;E" is therefore a valid value..

## We want to ignore cases involving B
dat <- subset(sheet1, !grepl("B", sheet1$OUTCOME))
## We also want to ignore cases which don't have an outcome
dat <- subset(dat, grepl("A|C", dat$OUTCOME))
dat$won <- as.numeric(grepl("A", dat$OUTCOME))


## ----getexperience-------------------------------------------------------
counsel.vars <- grep("APPCOUNS|RESPCOUNS", names(dat), value = TRUE)
experience.df <- melt(dat[,c("NEUTRAL", "CASEID","won", "HANDDOWN", counsel.vars)],
                      id.vars = c("NEUTRAL", "CASEID", "HANDDOWN", "won"))
experience.df <- subset(experience.df, !is.na(value))

normalize_names <- function(x) {
    x <- sub("\\s+QC", "", x)
    x <- sub("\\s+BL", "", x)
    x <- gsub("Dr |Prof |Sir ", "", x)
    x <- tolower(x)
    return(x)
    }

experience.df$normalForm <- normalize_names(experience.df$value)


### Remove duplicates within a single case
experience.df <- subset(experience.df,
                        !duplicated(experience.df[,c("NEUTRAL","normalForm")]))

### Create a running count of the experience
experience.df <- experience.df %>%
    group_by(normalForm) %>%
    arrange(HANDDOWN) %>%
    mutate(nAppearances = lag(1:(n())),
           nWins = lag(cumsum(won)))

experience.df$nAppearances[is.na(experience.df$nAppearances)] <- 0
experience.df$nWins[is.na(experience.df$nWins)] <- 0

experience.df$success <- (1 + experience.df$nWins) / (2 + experience.df$nAppearances)
experience.df$Side <- sub("\\d+", "", experience.df$variable)
### Aggregate by either side
experience.df <- experience.df %>%
    group_by(NEUTRAL, CASEID, Side) %>%
    arrange(variable) %>% 
    summarize(success = head(success, 1), # weighted.mean(success, (1 + nAppearances)),
              experience = head(nAppearances, 1)) # mean(nAppearances))

success.df <- dcast(experience.df, NEUTRAL + CASEID ~ Side,
                       value.var = "success")
names(success.df)[3:4] <- c("AppSuccess", "RespSuccess")
success.df$AppSuccess[is.na(success.df$AppSuccess)] <- 0.5
success.df$RespSuccess[is.na(success.df$RespSuccess)] <- 0.5

experience.df <- dcast(experience.df, NEUTRAL + CASEID ~ Side,
                       value.var = "experience")
names(experience.df)[3:4] <- c("AppExperience", "RespExperience")
experience.df$AppExperience[is.na(experience.df$AppExperience)] <- 0
experience.df$RespExperience[is.na(experience.df$RespExperience)] <- 0

dat <- merge(dat, experience.df, all = TRUE)
dat <- merge(dat, success.df, all = TRUE)


## ----win_plots-----------------------------------------------------------
terms <- c("Hil", "Eas", "Tri", "Mic")
factor.levels <- as.vector(t(outer(2009:2017, terms, paste, sep = " ")))

## ### Do by appellant type
type.recode.str <- "c('A','B')='Individual or association';
'C'='Company';
c('D','J')='Local govt., other public bodies';
c('E','F','H')='Central govt.';
else ='Unclear'"

dat$AppType.recoded <- car:::recode(dat$APPTYPE, type.recode.str)
dat$RespType.recoded <- car:::recode(dat$RESPTYPE, type.recode.str)

plot.df <- dat %>%
    mutate(year = as.numeric(format(HANDDOWN, "%Y")),
           term_short = substr(TERM, 0, 1),
           year_term = paste(year,
                             term_short,
                             sep = " ")) %>%
    mutate(year_term = factor(year_term,
                              levels = factor.levels,
                              ordered = TRUE)) %>%
    group_by(year, term_short, year_term) %>%
    summarize(date = as.Date(min(HANDDOWN)),
              WinRate = mean(won, na.rm = TRUE),
              nCases = sum(!is.na(won)))

plot.df$term_short <- factor(plot.df$term_short,
                             levels = c("M", "H", "E", "T"),
                             ordered = TRUE)

p1 <- ggplot(plot.df, aes(x = date, y = WinRate, weight = nCases)) +
    ## geom_smooth(se = FALSE, method = "lm", alpha = 0.5) + 
    geom_line(alpha = 0.3) +
    geom_point(color = "white", size = 6) + 
    geom_text(aes(label = term_short),
              size = 3) + 
    scale_y_continuous("Appellant success rate",
                       labels = scales::percent,
                       limits = c(0, 1)) +
    scale_x_date("Year and legal term") +
    scale_size_continuous(guide = FALSE) +
    ## geom_hline(yintercept = 0.5, col = "darkgrey") + 
    ## geom_point(aes(colour = term_short)) +
    theme_uksc()

p1.df <- plot.df

dat <- dat %>%
    group_by(Area) %>%
    mutate(nCases = n())

dat$Area_pretty <- paste0(dat$Area, "\n(", dat$nCases, ")")
dat$Area_pretty <- factor(dat$Area_pretty,
                          levels = unique(dat$Area_pretty),
                          ordered = TRUE)
dat$Area_pretty <- reorder(dat$Area_pretty,
                           dat$won)

p2 <- ggplot(subset(dat, !is.na(Area)), aes(x = Area_pretty, y = won, group = Area)) +
    stat_summary(fun.y = "mean", geom = "bar") +
    scale_x_discrete("Area of law") +
    scale_y_continuous("Appellant success rate", labels = scales::percent,
                       limits = c(0, 1)) +
    theme_uksc()

plot.df <- subset(dat, !is.na(OpinionBelow)) %>%
    group_by(OpinionBelow) %>%
    summarize(WinRate = mean(won),
              nCases = sum(!is.na(won)))

p3 <- ggplot(plot.df, aes(x = OpinionBelow, y = WinRate, size = nCases, weight = nCases)) +
    geom_point() +
    geom_smooth(se = FALSE, method = "lm") + 
    scale_x_continuous("Opinion below\n(higher values more favourable for appellant)",
                       limits = c(0.1, 0.5)) +
    scale_y_continuous("Appellant success rate",
                       labels = scales::percent,
                       limits = c(0, 1)) +
    scale_size_continuous(guide = FALSE) + 
    theme_uksc()

p3.df <- plot.df

plot.df <- subset(dat, AppType.recoded != "Unclear") %>%
    group_by(AppType.recoded) %>%
    summarize(WinRate = mean(won, na.rm = TRUE),
              nCases = sum(!is.na(won)))

dat <- dat %>%
    group_by(AppType.recoded) %>%
    mutate(nCases = n())

dat$AppType.pretty <- paste0(dat$AppType.recoded, "\n(", dat$nCases, " cases)") 
dat$AppType.pretty <- unlist(lapply(sapply(dat$AppType.pretty, function(x) {
    strwrap(x, width = 24)
    }), paste0, collapse = "\n"))

dat$AppType.pretty <- factor(dat$AppType.pretty,
                             levels = unique(dat$AppType.pretty),
                             ordered = TRUE)
dat$AppType.pretty <- reorder(dat$AppType.pretty,
                              dat$won)
                            
p4 <- ggplot(subset(dat, AppType.recoded != "Unclear"),
             aes(x = AppType.pretty, y = won, group = AppType.pretty)) +
    stat_summary(fun.y = "mean", geom = "bar") +
    scale_x_discrete("Appellant type") +
    scale_y_continuous("Appellant success rate", labels = scales::percent,
                       limits = c(0, 1)) +
    theme_uksc()



## ----p1, fig = TRUE, fig.cap = "Appellant success over time. Each point represents an average for a legal term (M = Michaelmas; H = Hilary; E = Easter; T = Trinity. ", fig.width = 6, fig.height = 4----
print(p1)


## ----p2, fig = TRUE, fig.cap = "Appellant success by area of law", fig.height = 4, fig.width = 6----
print(p2)


## ----p3, fig = TRUE, fig.cap = "Appellant success by opinion below. Each point represents an average for distinct values of the opinion below variable. Plotted points are scaled to represent the number of cases with that value. The solid blue line indicates the trend. ", fig.height = 4, fig.width = 6----
print(p3)


## ----p4, fig = TRUE, fig.cap = "Appellant success by type of appellant", fig.height = 4, fig.width = 6----
print(p4)


## ----mungech9------------------------------------------------------------
## ### Do by relative number of lawyers
applawyerpos <- grep("APPCOUNS", names(dat))
dat$nAppLawyers <- apply(dat[,applawyerpos], 1, function(x) sum(!is.na(x)))
dat$nAppLawyers[is.na(dat$nAppLawyers)] <- 0
resplawyerpos <- grep("RESPCOUNS", names(dat))
dat$nRespLawyers <- apply(dat[,resplawyerpos], 1, function(x) sum(!is.na(x)))
dat$nRespLawyers[is.na(dat$nRespLawyers)] <- 0
dat$relativeNoLawyers <- dat$nAppLawyers - dat$nRespLawyers


countQC <- function(x) sum(grepl("QC", x, ignore.case = TRUE))
dat$nAppQCs <- apply(dat[,applawyerpos], 1, countQC)
dat$nAppQCs[is.na(dat$nAppQCs)] <- 0
dat$nRespQCs <- apply(dat[,resplawyerpos], 1, countQC)
dat$relativeNoQCs <- dat$nAppQCs - dat$nRespQCs
dat$nRespQCs[is.na(dat$nRespQCs)] <- 0

dat$relativeNoIntervenors <- as.numeric(dat$APPINTER == "Yes") -
    as.numeric(dat$RESPINTER == "Yes")



dat$daysHearing <- apply(dat[,grep("HEARINGDATE", names(dat))],
                         1,
                         function(x)sum(!is.na(x)))
## Litigant types
## -	A. Natural person
## -	B. Association
## -	C. Company
## -	D. Local govt.
## -	E. National govt.
## -	F. Prosecuting or judicial authority
## -	G. Local govt. attorney-general, or similar
## -	H. National govt. attorney-general, or similar
## -	J. Other public body (school, hospital, police force) not part of nat'l govt

doStatus <- function(app, resp, codes) {
     retval <- rep(0, length(app))
     case1 <- which(is.element(app, codes) & !is.element(resp, codes))
     case2 <- which(!is.element(app, codes) & is.element(resp, codes))
     retval[case1] <- 1
     retval[case2] <- -1
     return(retval)
}
dat$isNP_Assoc <- doStatus(dat$APPTYPE, dat$RESPTYPE, c("A","B"))
dat$isCompany <- doStatus(dat$APPTYPE, dat$RESPTYPE, c("C"))
dat$isLocalGovt <- doStatus(dat$APPTYPE, dat$RESPTYPE, c("D", "G", "J"))
dat$isCentralGovt <- doStatus(dat$APPTYPE, dat$RESPTYPE, c("E", "H"))
dat$isJudAuth <- doStatus(dat$APPTYPE, dat$RESPTYPE, c("F"))

dat$CentralGovtIntervenes <- doStatus(dat$APPINTERTYPE,
                                      dat$RESPINTERTYPE,
                                      c("E", "H"))

### Did it go through PTA?
pta <- read.csv("./data/pta_data.csv")
pta$Date <- as.Date(pta$Date, "%d/%m/%Y")

### Make sure only to get decisions with a UKSC caseid
pta <- subset(pta, grepl("UKSC", pta$CASEID))

dat$wasPTA <- as.numeric(dat$CASEID %in% unique(pta$CASEID))

dat$Leapfrog <- as.numeric(is.na(dat$FICOURT))

dat$legalAided <- 0
dat$legalAided[which(dat$APPASSIST == "Yes" & dat$RESPASSIST == "No")] <- 1
dat$legalAided[which(dat$APPASSIST == "No" & dat$RESPASSIST == "Yes")] <- -1


## ----quickintervenormunge------------------------------------------------
dat$relativeNoIntervenors[is.na(dat$relativeNoIntervenors)] <- 0


## ----finalmunge----------------------------------------------------------
dat$Area <- relevel(dat$Area, "Civil")
actor.vars <- c("isNP_Assoc","isCompany","isLocalGovt","isCentralGovt",
          "isJudAuth")
### At some point, will need to introduce changes for Scottish appeals
### This is due to changes made in 2014/5
### See https://www.supremecourt.uk/news/new-powers-for-Scotlands-appeal-court-in-force-today.html
dat$postCutoff <- as.numeric(dat$APPORDDATE > as.Date("2015-09-22"))
dat$postCutoff[is.na(dat$postCutoff)] <- 0
dat$postCutoffScots <- dat$postCutoff * as.numeric(dat$Area == "Scots")
dat$OpinionBelow[is.na(dat$OpinionBelow)] <- 0.5


dat$pub_or_chan <- as.numeric(dat$Area %in% c("Public", "Chancery"))


vars <- c("OpinionBelow",
          "Importance",
          "relativeNoLawyers",
          "relativeNoQCs",
          "relativeNoIntervenors")
            
means_std_vars <- numeric()
sd_std_vars <- numeric()
for (v in vars) {
    means_std_vars <- c(means_std_vars,
                        mean(dat[[v]], na.rm = TRUE))
    sd_std_vars <- c(sd_std_vars,
                        sd(dat[[v]], na.rm = TRUE))
    dat[,v] <- arm::rescale(dat[[v]])
}

names(means_std_vars) <- names(sd_std_vars) <- vars

keep.vars <- c("won",
               "OpinionBelow",
               "Area",
               "postCutoff",
               "postCutoffScots",
               "wasPTA",
               "relativeNoLawyers",
               "relativeNoQCs",
               "relativeNoIntervenors",
               "legalAided",
               "isNP_Assoc",
               "isLocalGovt",
               "isCentralGovt",
               "CentralGovtIntervenes",
               "pub_or_chan")

dat$Area <- relevel(dat$Area,
                    "Civil")
dat$Area[is.na(dat$Area)] <- "Civil"
dat$postCutoffScots[is.na(dat$postCutoffScots)] <- 0



## ----modellingch9--------------------------------------------------------
if (file.exists("cache/09_cached_results.RData")) {
    load("cache/09_cached_results.RData")
} else { 
    legalmod <- stan_glm(won ~ OpinionBelow +
                             Area +
                             postCutoff + 
                             postCutoffScots +
                             wasPTA,
                         family = binomial,
                         data = dat,
                         cores = 2, chains = 2, seed = 1982)
    
    resourcemod <- stan_glm(won ~ relativeNoLawyers +
                                relativeNoQCs +
                                relativeNoIntervenors +
                                legalAided,
                            family = binomial,
                            data = dat,
                            cores = 2, chains = 2, seed = 1982)
    
    politicsmod <- stan_glm(won ~ isNP_Assoc +
                                isLocalGovt +
                                (isCentralGovt +
                                 CentralGovtIntervenes) * 
                                pub_or_chan,
                            family = binomial,
                            data = dat,
                            cores = 2, chains = 2, seed = 1982)
    
    combmod <- stan_glm(won ~ OpinionBelow + Area +
                            postCutoff +
                            postCutoffScots + 
                            wasPTA +
                            relativeNoLawyers + relativeNoQCs +
                            relativeNoIntervenors + legalAided +
                            isNP_Assoc + isLocalGovt +
                            isCentralGovt + CentralGovtIntervenes +
                            isCentralGovt:pub_or_chan + CentralGovtIntervenes:pub_or_chan,
                        data = dat,
                        family = binomial,
                        cores = 2, chains = 2, seed = 1982)
    save(combmod, politicsmod, resourcemod, legalmod,
         file = "cache/09_cached_results.RData")
    }



## ----legalmodch9, fig = TRUE, fig.cap = "Legal model of appellate outcomes. Thick lines indicate 83 percent credible intervals. Thin lines indicate 95 percent credible intervals. ", fig.width = 6, fig.height = 4----
legalmod.labs <- list(
'(Intercept)'= 'Intercept',
'scale(OpinionBelow)'='Opinion below',
'AreaPublic'= 'Public law case',
'AreaChancery' = 'Tax/Chancery case',
'AreaCriminal' = 'Criminal case',
'AreaFamily' = 'Family law case',
'AreaN Ireland' = 'Northern Irish case',
'AreaScots' = 'Scots case',
'postCutoff' = 'Case heard after Sep. 2015',
'wasPTA' = 'Case received PTA',
'postCutoffScots' = 'Scots case heard after Sep. 2015') 


mycoefplot(legalmod,
           ylab = "Coefficient value",
           variable.labels = legalmod.labs,
           sec.ylab = "Odds are ... times greater")

family.coef <- round(coef(legalmod)["AreaFamily"], 3)
family.fx <- round(exp(coef(legalmod)["AreaFamily"]), 3)

opbelow.coef <- round(coef(legalmod)["OpinionBelow"], 3)
opbelow.fx <- round(exp(coef(legalmod)["OpinionBelow"] * 1/6), 3)

legal.pcp <- hitmiss(legalmod)
comb.pcp <- hitmiss(combmod)


## ----orgmodoutch9, fig = TRUE, fig.cap = "Resource-based model of appellate outcomes. Thick lines indicate 83 percent credible intervals. Thin lines indicate 95 percent credible intervals. ", fig.width = 6, fig.height = 2.5----
orgmod.labs <- list(
'(Intercept)'='Intercept',
'relativeNoLawyers'='Relative no. of lawyers', 
'relativeNoQCs'='Relative no. of QCs',
'relativeNoIntervenors'='Relative no. of interveners',
'legalAided'='Appellant legally aided')

mycoefplot(resourcemod,
           ylab = "Coefficient value",
           variable.labels = orgmod.labs,
           sec.ylab = "Odds are ... times greater")



## ----polmodoutch9, fig = TRUE, fig.cap = "Political model of appellate outcomes. Thick lines indicate 83 percent credible intervals. Thin lines indicate 95 percent credible intervals. ", fig.width = 6, fig.height = 3.5----
polmod.labs <- list(
'(Intercept)'='Intercept', 
'isNP_Assoc'='Natural person/assn.',
'isLocalGovt'='Local govt',
'isCentralGovt'='Central govt',
'CentralGovtIntervenes'='Central govt (intervenor)',
'pub_or_chan'='Public law/tax case',
'isCentralGovt:pub_or_chan'='Central govt × public/tax case',
'CentralGovtIntervenes:pub_or_chan'='Central govt (intervenor) × public/tax case')

mycoefplot(politicsmod,
           ylab = "Coefficient value",
           variable.labels = polmod.labs,
           sec.ylab = "Odds are ... times greater")



## ----combmodoutch9, fig = TRUE, fig.cap = "Combined model of appellate outcomes. Thick lines indicate 83 percent credible intervals. Thin lines indicate 95 percent credible intervals. ", fig.width = 6, fig.height = 8----
combmod.labs <- c(legalmod.labs, orgmod.labs, polmod.labs)
combmod.labs <- combmod.labs[!duplicated(combmod.labs)]
### Add on interaction terms with reverse order
addons <- combmod.labs
names(addons) <- gsub("(.*):(.*)", "\\2:\\1", names(combmod.labs))
combmod.labs <- c(combmod.labs, addons)

mycoefplot(combmod,
           variable.labels = combmod.labs,
           sec.ylab = "Odds are ... times greater",
           ylab = "Coefficient value")





## ----discretechangesch9--------------------------------------------------
predictprob <- function(mod, newdata = NULL) {
    pr <- posterior_linpred(mod, newdata = newdata,
                            transform = TRUE)
    pr[,1]
}

cf <- function(mod, old, new, alpha = 0.05) {
    oldprobs <- predictprob(mod, newdata = old)
    newprobs <- predictprob(mod, newdata = new)
    oldprobs <- oldprobs
    newprobs <- newprobs
    diff <- newprobs - oldprobs
    return(data.frame(mean = median(diff),
                      llo = quantile(diff, 1/12),
                      hhi = quantile(diff, 11/12),
                      lo = quantile(diff, alpha/2),
                      hi = quantile(diff, 1 - alpha / 2)))
}

civil_case <- with(dat,
                   data.frame(OpinionBelow = mean(OpinionBelow, na.rm = TRUE),
                              Area = "Civil",
                              postCutoff = 0,
                              postCutoffScots = 0,
                              wasPTA = 0,
                              relativeNoLawyers = 0,
                              relativeNoQCs = 0,
                              relativeNoLawyers = 0,
                              relativeNoIntervenors = 0,
                              legalAided = 0,
                              isNP_Assoc = 0,
                              isLocalGovt = 0,
                              isCentralGovt = 0,
                              CentralGovtIntervenes = 0,
                              pub_or_chan = 0))

public_case <- civil_case
public_case$Area <- "Public"
public_case$pub_or_chan <- 1

### a) Opinion below
inc_opinion_civil <- civil_case
inc_opinion_public <- public_case
inc_opinion_civil$OpinionBelow <- inc_opinion_civil$OpinionBelow +
    2 * sd(dat$OpinionBelow, na.rm = TRUE)
inc_opinion_public$OpinionBelow <- inc_opinion_public$OpinionBelow +
    2 * sd(dat$OpinionBelow, na.rm = TRUE)

opinion.chg.civil <- cf(combmod, civil_case, inc_opinion_civil)
opinion.chg.public <- cf(combmod, public_case, inc_opinion_public)

### b) Relative no. QCs
inc_qcs_civil <- civil_case
inc_qcs_public <- public_case
inc_qcs_civil$relativeNoQCs <- 1
inc_qcs_public$relativeNoQCs <- 1

qc.chg.civil <- cf(combmod, civil_case, inc_qcs_civil)
qc.chg.public <- cf(combmod, public_case, inc_qcs_public)

### c) PTA
inc_pta_civil <- civil_case
inc_pta_public <- public_case
inc_pta_civil$wasPTA <- 1
inc_pta_public$wasPTA <- 1

pta.chg.civil <- cf(combmod, civil_case, inc_pta_civil)
pta.chg.public <- cf(combmod, public_case, inc_pta_public)

### (d) Central govt.
inc_govt_civil <- civil_case
inc_govt_public <- public_case

inc_govt_civil$isCentralGovt <- 1
inc_govt_public$isCentralGovt <- 1

govt.chg.civil <- cf(combmod, civil_case, inc_govt_civil)
govt.chg.public <- cf(combmod, public_case, inc_govt_public)

out <- rbind(opinion.chg.civil, #opinion.chg.public,
            qc.chg.civil, #qc.chg.public,
             pta.chg.civil, #pta.chg.public,
             govt.chg.civil, govt.chg.public)

out <- as.data.frame(out)
changes <- c("3/3 judges rule against appellant to 3/4 judges rule against appellant",
"One additional QC for appellant",
"Case received PTA",
"Company to central govt. appellant",
"Company to central govt. appellant, public law case")

changes <- str_wrap(changes, 20)
out$Change <- changes



## ----predchangesch9, fig = TRUE, fig.cap = "Changes in the predicted probability of appellant success, according to certain specific changes. The baseline is a company appealing in a civil case which did not go through the permission to appeal process. ", fig.width = 4, fig.height = 6----

dodge.width <- 0.9

## Change bar colors
out$isSignificant <- apply(out[,c("lo", "hi")], 1, isSig)

out$inner.bar.col <- ifelse(out$isSignificant,
                            inner.bar.col,
                            "lightgrey")
out$outer.bar.col <- ifelse(out$isSignificant,
                            outer.bar.col,
                            "darkgrey")

ggplot(out, aes(x = Change,
                y = mean, ymin = lo, ymax = hi)) +
    scale_x_discrete("") +
    scale_y_continuous("Change in probability\nofappellant success\n(percentage points)",
                       labels = scales::percent) +
    geom_linerange(aes(colour = inner.bar.col),
                   position = position_dodge(width = dodge.width),
                   size=inner.bar.size, 
                   alpha = 0.8) +
    geom_linerange(aes(ymin = llo, ymax = hhi,
                       colour = outer.bar.col),
                   position = position_dodge(width = dodge.width),
                   size = outer.bar.size,
                   alpha = 0.8) +
    geom_hline(aes(x = 0, yintercept = 0), lty = 1) +
    geom_point(aes(colour = outer.bar.col),
               shape = 1,
               position = position_dodge(width = dodge.width),
               size= (inner.bar.size + outer.bar.size) * 2.5,
               stroke = 1.5) +
    scale_colour_identity() +
    scale_fill_identity() +
    theme_uksc() +
    coord_flip() +
    theme(legend.position = "bottom",
          legend.dir = "horizontal")




## ----fitstatsch9, fig = TRUE, fig.cap = "Goodness-of-fit statistics for different models. Thin lines indicate 95 percent credible intervals; thick lines indicate 83 percent credible intervals. ", warning = FALSE----

epcp.rstanarm <- function(mod) {
    
    poslp <- posterior_linpred(mod)
    sumdf <- as.data.frame(summary(mod)) 
    preds <- (poslp > 0.5)
    
    actual <- matrix(mod$y,
                     nrow = nrow(poslp),
                     ncol = length(mod$y),
                     byrow = TRUE)
    correct <- (preds == actual)

    correct <- apply(correct, 1, mean)
    epcp <- data.frame(term = "PCP\n(bigger is better)",
                       mean = mean(correct),
                       lo = quantile(correct, 0.025),
                       hi = quantile(correct, 0.975),
                       llo = quantile(correct, 1/12),
                       hhi = quantile(correct, 11/12))

    the.loo <- loo(mod)
    the.loo <- data.frame(term = "LOOIC\n(smaller is better)",
                          mean = the.loo$looic,
                          lo = the.loo$looic + qnorm(1/40) * the.loo$se_looic,
                          hi = the.loo$looic + qnorm(39/40) * the.loo$se_looic,
                          llo = the.loo$looic + qnorm(1/24) * the.loo$se_looic,
                          hhi = the.loo$looic + qnorm(23/24) * the.loo$se_looic)
                                        #    return(rbind(epcp, the.loo))
    return(rbind(epcp, the.loo))
}

if (file.exists("cache/09_gof.RData")) {
    load("cache/09_gof.RData")
} else { 
    

    legal.stats <- epcp.rstanarm(legalmod)
    org.stats <- epcp.rstanarm(resourcemod)
    pol.stats <- epcp.rstanarm(politicsmod)
    comb.stats <- epcp.rstanarm(combmod)

    legal.stats$Model <- "Legal"
    org.stats$Model <- "Organisational"
    pol.stats$Model <- "Political"
    comb.stats$Model <- "Combined"

    null.stats <- data.frame(Model = "Null",
                             term = "PCP\n(bigger is better)",
                             ## mean = mean(combmod.rstanarm$y == 1, na.rm = TRUE))
                             mean = mean(legalmod$y == 1, na.rm = TRUE))
    
    gof.df <- rbind(
        rbind(
            rbind(legal.stats,
                  org.stats),
            pol.stats),
        comb.stats)
    
    gof.df <- merge(gof.df, null.stats, all = TRUE)
    
    gof.df$Model <- factor(gof.df$Model,
                           levels = rev(c("Null", "Legal",
                                          "Organisational", "Political",
                                          "Combined")),
                           ordered = TRUE)

    legal.loo <- loo(legalmod)
    combined.loo <- loo(combmod)
    political.loo <- loo(politicsmod)
org.loo <- loo(resourcemod)

    save(list = c("legal.loo", "combined.loo",
                  "political.loo", "org.loo",
"gof.df"),
         file = "cache/09_gof.RData")


}

ggplot(gof.df, aes(x = Model, y = mean, ymax = hi, ymin = lo)) +
    geom_linerange(size = inner.bar.size,
                   colour = inner.bar.col,
                   alpha = 0.8) +
        geom_linerange(aes(ymin = llo, ymax = hhi),
                       size = outer.bar.size,
                       colour = outer.bar.col,
                       alpha = 0.8) +
    geom_point(aes(colour = outer.bar.col),
               size= (inner.bar.size + outer.bar.size) * 2.5,
               shape= my.pchs[1],
               stroke = 1.5) +
    scale_colour_identity() + 
    facet_wrap(~term, scales = "free_x") +
    scale_y_continuous("Value") +
    coord_flip() +
    theme_uksc()

combined.pcp <- subset(gof.df,
                       Model == "Combined" & grepl("PCP", term)) %>%
    pull(mean)

combined.pcp <- round(100 * combined.pcp)

